library(stargazer)
library(survival)
library(MASS)
library(tidyverse)
library(ggplot2)
library(reshape2)
library(ggthemes)
library(CFC)
library(SurvRegCensCov)
library(mlogit)
library(sandwich)
library(lmtest)
library(BMA)

source("utility_functions_WRVC.R")
source("BMA_function_aicsurv_adjusted.R")

dat <- read_csv("Who_Restarts_Violent_Conflict_Replication.csv")

main_dat <- arrangement(dat, preact_def = 2, recast_def = 1, both = 1)

##################################
#### Figure 1 in the main text####
##################################

decades <- c(1960, 1970, 1980, 1990, 2000, 2010)
decades_out <- data.frame(decade = decades, repeated = NA, recast = NA)
for (i in decades) {
  p <- subset(main_dat, main_dat$year >= i & main_dat$year < i + 10)
  decades_out[decades_out$decade == i, ]$repeated <- sum(p$status == 1)
  decades_out[decades_out$decade == i, ]$recast <- sum(p$status == 2)
}

for_ts_plot <- melt(decades_out, id.vars = "decade")

ts_g <- ggplot(for_ts_plot, aes(x = decade, y = value, fill = variable)) + theme_classic() +
  geom_bar(position = "dodge", stat = "identity") + labs(fill = "types") + ylab("Count") + 
  scale_x_continuous(breaks = c(1960, 1970, 1980, 1990, 2000, 2010),
                     labels = c("1960s", "1970s", "1980s", "1990s", "2000s", "2010s")) +
  scale_fill_manual(values = c("#66CCFF", "#CC0033"))
ts_g
## To save
#ggsave("time_series_reprec.png", width = 6, height = 4)

#################################################################
#### Main Analysis (Table 3 in the main text & Appendix E.6) ####
#################################################################

Full_Model <- "Num_Actor + Revolution + Res_Conf + PKO + Mountain + lag_IMR + ELF + Victory + lag_Growth + Agreement + Post_CW"

main_res <- competing_risks(main_dat, indep = Full_Model, subest = T)

## Extracting dataframe used for the main analysis. This will be used later for descriptive statistics.
compdat <- model_extract(main_dat, Full_Model)

## Table 3 in the main text
main_table <- stargazer(main_res[[1]], main_res[[4]], main_res[[5]], style = "aer",
                  omit.stat = c("rsq", "max.rsq", "wald", "logrank"),
                  font.size = "footnotesize", digits = 2,
                  covariate.labels = c("Number of Groups", "Revolution",
                               "Resource Conflict", "PKO",
                               "Mountain", "lag IMR", "ELF", "Victory", "lag Growth",
                               "Peace Agreement", "Post Cold War"))
## To save
#cat(main_table, sep = "\n", file = "main_table.tex")

## The numbers of events was added manually
freq_tab(compdat)

## Results from subdistribution models (Appendix E.6)
### Table 11
main_subdis_table <- stargazer(main_res[[1]], main_res[[2]], main_res[[3]], style = "aer",
                        omit.stat = c("rsq", "max.rsq", "wald", "logrank"),
                        font.size = "footnotesize", digits = 2,
                        covariate.labels = c("Number of Groups", "Revolution",
                                             "Resource Conflict", "PKO",
                                             "Mountain", "lag IMR", "ELF", "Victory", "lag Growth",
                                             "Peace Agreement", "Post Cold War"))
## To save
#cat(main_subdis_table, sep = "\n", file = "main_subdis_table.tex")


#############################################
#### Descriptive Statistics (Appendix C) ####
#############################################

unit_data <- subset(compdat, compdat$t_time==1)
descri <- unit_data[,c(8:12, 15, 17, 26, 31, 36, 38, 43, 45, 52, 54, 55, 57, 58, 60, 63, 64)]
stargazer(descri)

### KM plots

### Pooled KM plot (Figure 3)
png("KMP_P.png", width = 400, height = 300)
plot(survfit(Surv(compdat$t_time-1, compdat$t_time, compdat$event) ~ 1, data = compdat),
     xlab = "time", ylab = "Survival Function")
dev.off()

### KM plot for repeated conflict (Figure 4)
png("KMP_Rep.png", width = 400, height = 300)
plot(survfit(Surv(compdat$t_time-1, compdat$t_time, compdat$status == 1) ~ 1, data = compdat),
     xlab = "time", ylab = "Survival Function")
dev.off()

### KM plot for recast conflict (Figure 5)
png("KMP_Rec.png", width = 400, height = 300)
plot(survfit(Surv(compdat$t_time-1, compdat$t_time, compdat$status == 2) ~ 1, data = compdat),
     xlab = "time", ylab = "Survival Function")
dev.off()

#########################################################################################
#### Excluding irrelevant conflict from recurrence using Zeigler 2016 (Appendix E.1) ####
#########################################################################################

z1_dat <- subset(dat, dat$Zeigler_Exclude == 0)
z2_dat <- subset(dat, dat$Zeigler_Extended == 0)

z1 <- arrangement(z1_dat, preact_def = 2)
z2 <- arrangement(z2_dat, preact_def = 2)

z1_res <- competing_risks(z1, Full_Model)
z2_res <- competing_risks(z2, Full_Model)

### Table 6
zeigler_table <- stargazer(z1_res[[1]], z1_res[[4]], z1_res[[5]], z2_res[[1]], z2_res[[4]], z2_res[[5]],
                           style = "aer",
                           omit.stat = c("rsq", "max.rsq", "wald", "logrank"),
                           font.size = "footnotesize", digits = 2,
                           covariate.labels = c("Number of Groups", "Revolution",
                               "Resource Conflict", "PKO",
                               "Mountain", "lag IMR", "ELF", "Victory", "lag Growth",
                               "Peace Agreement", "Post Cold War"))
## To save
#cat(zeigler_table, sep = "\n", file = "zeigler_table.tex")

#################################################################
#### Splinter Recurrence as Repeated Conflict (Appendix E.2) ####
#################################################################

### Figure 17

spl_plot_dat <- subset(main_dat, main_dat$status > 0)
barp_dat <- data.frame(table(spl_plot_dat$state, spl_plot_dat$Spl_Rec))
names(barp_dat) <- c("RecurType", "BySplinter", "Frequency")
barp_dat$RecurType <- sapply(barp_dat$RecurType, function(x){ifelse(x == "1", "Repeat", "Recast")})
barp_dat$BySplinter <- sapply(barp_dat$BySplinter, function(x){ifelse(x == "0", "No", "Yes")})

ggplot(barp_dat, mapping = aes(x = RecurType, y = Frequency, fill = BySplinter)) + theme_classic() +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("#66CCFF", "#CC0033"))
## To save
#ggsave("Splinter_Recurrence.jpg", width = 6, height = 4)

### Competing risk analysis coding recurrence by splinters as repeated conflict

spl_rep_dat <- main_dat
re_status <- spl_rep_dat$status
for (i in 1:nrow(spl_rep_dat)) {
  if(re_status[i] == 2 & spl_rep_dat$Spl_Rec[i] == 1){
    re_status[i] <- 1
  }
}
spl_rep_dat$status <- re_status

splinter_res <- competing_risks(spl_rep_dat, indep = Full_Model, subest = F)

### Table 7
spl_table <- stargazer(splinter_res[[1]], splinter_res[[4]], splinter_res[[5]], style = "aer",
                       omit.stat = c("rsq", "max.rsq", "wald", "logrank"),
                       font.size = "footnotesize", digits = 2,
                       covariate.labels = c("Number of Groups", "Revolution",
                                            "Resource Conflict", "PKO",
                                            "Mountain", "lag IMR", "ELF", "Victory", "lag Growth",
                                            "Peace Agreement", "Post Cold War"))
## To save
#cat(spl_table, sep = "\n", file = "spl_table.tex")

######################################################
#### Rebel victory cases included  (Appendix E.3) ####
######################################################

rv_dat <- arrangement(dat, preact_def = 2, RV = 1)

rv_res <- competing_risks(rv_dat, Full_Model, subest = F)

### Table 8
rv_table <- stargazer(rv_res[[1]], rv_res[[4]], rv_res[[5]], style = "aer",
                       omit.stat = c("rsq", "max.rsq", "wald", "logrank"),
                       font.size = "footnotesize", digits = 2,
                       covariate.labels = c("Number of Groups", "Revolution",
                                            "Resource Conflict", "PKO",
                                            "Mountain", "lag IMR", "ELF", "Victory", "lag Growth",
                                            "Peace Agreement", "Post Cold War"))
## To save
#cat(rv_table, sep = "\n", file = "rv_table.tex")

##########################################
##### Other Intervals (Appendix E.4) #####
##########################################

oi_dat <- main_dat %>% group_by(pid) %>% mutate(max_ttime = max(t_time))

oi_spl <- split(oi_dat, oi_dat$cid)

## Erase peace with duration less than 3 years
short_peace_erase <- function(x, interval){
  this <- x %>% group_by(pid) %>% mutate(lastyear = max(year))
  
  for(id in unique(this$pid)){
    if(this[this$pid == id, "max_ttime"][1,1] < interval){
      res_dur <- this[this$pid == id, "max_ttime"][1,1] + 1
      ly_omit <- this[this$pid == id, "lastyear"][1,1]
      this <- this[this$pid != id,]
      
      if(nrow(this) > 0){
        ## log Duration variable is adjusted accordingly
        if(sum(this$lastyear > ly_omit) > 0){
          target <- this$lastyear == min(this$lastyear[this$lastyear > ly_omit])
          this[target, "ln_Duration"] <- log(exp(unlist(c(this[target, "ln_Duration"]))) + res_dur)
        }
      }
    }
  }
  return(this)
}

### 2 years

oi2_spl <- lapply(oi_spl, short_peace_erase, interval = 3)

oi2 <- oi2_spl[[1]]
for (i in 2:length(oi2_spl)){
  oi2 <- rbind(oi2, oi2_spl[[i]])
}

oi2_res <- competing_risks(oi2, indep = Full_Model, subest = F)

compoi2 <- model_extract(oi2, Full_Model)
freq_tab(compoi2)

### 3 years

oi3_spl <- lapply(oi_spl, short_peace_erase, interval = 4)

oi3 <- oi3_spl[[1]]
for (i in 2:length(oi3_spl)){
  oi3 <- rbind(oi3, oi3_spl[[i]])
}

oi3_res <- competing_risks(oi3, indep = Full_Model, subest = F)

compoi3 <- model_extract(oi3, Full_Model)
freq_tab(compoi3)

### Table 9
int23_table <- stargazer(oi2_res[1], oi2_res[4], oi2_res[5], oi3_res[[1]], oi3_res[[4]], oi3_res[[5]], style = "aer",
                         omit.stat = c("rsq", "max.rsq", "wald", "logrank"),
                         font.size = "footnotesize", digits = 2,
                         covariate.labels = c("Number of Groups", "Revolution",
                                              "Resource Conflict", "PKO",
                                              "Mountain", "lag IMR", "ELF", "Victory", "lag Growth",
                                              "Peace Agreement", "Post Cold War"))
#cat(int23_table, sep = "\n", file = "int23_table.tex")

#######################################
#### Big/Small Wars (Appendix E.5) ####
#######################################

big <- subset(main_dat, main_dat$cumintense == 1)

res_big <- competing_risks(big, indep = Full_Model, subest = F)

compdat_big <- model_extract(big, Full_Model)

small <- subset(main_dat, main_dat$cumintense == 0)

res_small <- competing_risks(small, indep = Full_Model, subest = F)

compdat_small <- model_extract(small, Full_Model)
length(unique(compdat_small$pid))

### Table 10
sbwar_table <- stargazer(res_big[[1]], res_big[[4]], res_big[[5]], res_small[[1]], res_small[[4]], res_small[[5]],
                         style = "aer",
                         omit.stat = c("rsq", "max.rsq", "wald", "logrank"),
                         font.size = "footnotesize", digits = 2,
                         covariate.labels = c("Number of Groups", "Revolution",
                                              "Resource Conflict", "PKO",
                                              "Mountain", "lag IMR", "ELF", "Victory", "lag Growth",
                                              "Peace Agreement", "Post Cold War"))
#cat(sbwar_table, sep = "\n", file = "sbwar_table.tex")
freq_tab(compdat_big)
freq_tab(compdat_small)

#####################################################
#### Competing risk Weibull model (Appendix E.5) ####
#####################################################

weibull_res <- cfc.survreg(formula(paste("Surv(t_time-1, t_time, status) ~", Full_Model)), main_dat)
weibull_p <- survreg(formula(paste("Surv(t_time, event) ~", Full_Model)), main_dat)

con_resp <- ConvertWeibull(weibull_p)$vars
con_res1 <- ConvertWeibull(weibull_res$regs[[1]])$vars
con_res2 <- ConvertWeibull(weibull_res$regs[[2]])$vars

rresp <- round(con_resp, 2)
rres1 <- round(con_res1, 2)
rres2 <- round(con_res2, 2)

### Table 12
w_table <- stargazer(weibull_p, weibull_res$regs[[1]], weibull_res$regs[[2]], style = "aer",
                     omit.stat = c("rsq", "max.rsq", "wald", "logrank"),
                     font.size = "footnotesize", digits = 2,
                     covariate.labels = c("Number of Groups", "Revolution",
                                          "Resource Conflict", "PKO",
                                          "Mountain", "lag IMR", "ELF", "Victory", "lag Growth",
                                          "Peace Agreement", "Post Cold War"),
                     add.lines = list(c("shape", paste0(rresp[2,1], "$^{*}$"), paste0(rres1[2,1], "$^{***}$"), paste0(rres2[2,1], "$^{*}$")),
                                      (c("", paste0("(", rresp[2,2], ")"), paste0("(", rres1[2,2], ")"), paste0("(", rres2[2,2], ")")))))
## To save
#cat(w_table, sep = "\n", file = "weibull_table.tex")

#########################################
#### Multinomial Logit (Appendix E5) ####
#########################################

mldat <- main_dat

mldat$t_time2 <- mldat$t_time^2
mldat$t_time3 <- mldat$t_time^3
mldat$cid <- as.character(mldat$cid)

mlogit_dat <- mlogit.data(mldat, choice = "status", shape = "wide")

mres <- mlogit(formula = formula(paste0("status ~ 0| ", Full_Model, "+ t_time + t_time2 + t_time3|0")), data = mlogit_dat)

## A function to calculate the clustered standard errors for multinomial logit
cl.mlogit <- function(fm, cluster){
  
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- length(coefficients(fm))
  # edit 6/22/2015: change dfc
  # dfc <- (M/(M-1))*((N-1)/(N-K))
  dfc <- (M/(M-1))
  uj <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat.=crossprod(uj)/N)
  coeftest(fm, vcovCL)
}

ml_model <- model_extract(mldat, Full_Model)
cluster <- cl.mlogit(fm = mres, cluster = ml_model$cid)

### Table 13
ml_table <- stargazer(cluster, style = "ajps", digits = 2)
## To save
#cat(ml_table, sep = "\n", file = "ml_table.tex")

###############################################
#### Bayesian Model Averaging (Appendix F) ####
###############################################

## A set of variables to be used for BMA
whole_vars <- c("Victory", "Agreement", "PKO", "Revolution", "ln_Deaths", "ln_Duration", "Num_Actor",
                "Res_Conf", "Fragmentation",
                "ELF", "Mountain", "Post_CW", "lag_ln_GDPpc", "lag_ln_Pop", "lag_Growth",
                "lag_Democ", "lag_Autoc", "lag_IMR", "lag_LDG", "lag_Conf_Same",
                "lag_Conf_Neighb")

naomit <- na.omit(main_dat[,c(whole_vars, "year", "cid", "pid", "event", "t_time", "country",
                        "status", "state")])

mini <- naomit[, c(whole_vars, "t_time", "status")]
bma1 <- bic.surv(formula(Surv(time = t_time, event = status == 1) ~ .), data = mini)
bma2 <- bic.surv(formula(Surv(time = t_time, event = status == 2) ~ .), data = mini)

bma_aic1 <- aic.surv(formula(Surv(time = t_time, event = status == 1) ~ .), data = mini)
bma_aic2 <- aic.surv(formula(Surv(time = t_time, event = status == 2) ~ .), data = mini)

## Put 1 here to plot and save the output
  
if(0){
### Figure 18-20
pdf("bma_aic1.pdf")
plot.bic.surv(bma_aic1, mfrow = c(4,2), ylim = c(0,1), ask = "N")
dev.off()

## Figure 21-23
pdf("bma_aic2.pdf")
plot.bic.surv(bma_aic2, mfrow = c(4,2), ylim = c(0,1))
dev.off()

## Figure 24-26
pdf("bma_bic1.pdf")
plot.bic.surv(bma1, mfrow = c(4,2), ylim = c(0,1))
dev.off()

## Figure 27-29
pdf("bma_bic2.pdf")
plot.bic.surv(bma2, mfrow = c(4,2), ylim = c(0,1))
dev.off()
}

#########################
#### End of the code ####
#########################

